home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / toolbox / ode78.r < prev    next >
Encoding:
Text File  |  1994-07-18  |  4.4 KB  |  160 lines  |  [TEXT/ttxt]
  1. // ODE78  Integrates a system of ordinary differential equations using
  2. // 7th order formulas.
  3. //
  4. // out = ode78(F, t0, tfinal, y0, Uoutput, tol, trace)
  5. //
  6. // INPUT:
  7. // Ufun  - variable containing user-supplied problem description.
  8. //         Call: yprime = fun(t,y)
  9. //         t      - Time (scalar).
  10. //         y      - Solution row-vector.
  11. //         yprime - Returned derivative vector; yprime(i) = dy(i)/dt.
  12. // t0    - Initial value of t.
  13. // tfinal- Final value of t.
  14. // y0    - Initial value column-vector.
  15. // Uoutput User output function
  16. // tol   - The desired accuracy. (Default: tol = 1.e-6).
  17. // trace - If nonzero, each step is printed. (Default: trace = 0).
  18. //
  19. // OUTPUT:
  20. // out  - Returned integration time points (1st column).
  21. //        Remaining columns are solutions of equations specified in Ufun
  22. //
  23. // The result can be displayed by: plot(tout, yout).
  24.  
  25. //   Daljeet Singh
  26. //   Dept. Of Electrical Engg., The University of Alabama.
  27. //   11-24-1988.
  28.  
  29. //   Modified by Ian Searle
  30. //   4-18-93
  31.  
  32. static (collapse);
  33.  
  34. ode78 = function (Ufun, t0, tfinal, y0, Uoutput, tol, trace)
  35. {
  36.   local (alpha, beta, chi, delta, f, gamma1, h, hmax, ...
  37.          hmin, j, psi, pow, t, tau, tout, y, yout, lout, I, output);
  38.  
  39.   // The Fehlberg coefficients:
  40.   alpha = [ 2./27., 1/9,  1/6,    5/12,      .5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
  41.   beta =  [ 2/27,    0,     0,       0,       0, 0, 0, 0, 0, 0, 0, 0, 0 ;
  42.       1/36, 1/12,     0,       0,       0, 0, 0, 0, 0, 0, 0, 0, 0 ;
  43.       1/24,    0,   1/8,       0,       0, 0, 0, 0, 0, 0, 0, 0, 0 ;
  44.       5/12,    0,-25/16,   25/16,       0, 0, 0, 0, 0, 0, 0, 0, 0 ;
  45.        .05,    0,     0,     .25,      .2, 0, 0, 0, 0, 0, 0, 0, 0 ;
  46.    -25/108,    0,     0, 125/108,  -65/27, 125/54,  0, 0, 0, 0, 0, 0, 0 ;
  47.     31/300,    0,     0,       0,  61/225,    -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
  48.          2,    0,     0,   -53/6,  704/45,  -107/9,  67/90, 3, 0, 0, 0, 0, 0 ;
  49.    -91/108,    0,     0,  23/108,-976/135,  311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
  50.  2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0,...
  51.                                                                        0, 0 ;
  52.   3/205, 0, 0,  0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
  53.   -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, ...
  54.   51/82, 33/164, 12/41, 0, 1, 0 ]';
  55.  
  56.   chi = [0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';
  57.   psi = [1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1, -1,  -1 ]';
  58.   pow = 1/8;
  59.   if (!exist (tol)) { tol = 1.e-6; }
  60.   if (!exist (trace)) { trace = 0; }
  61.  
  62.   // Initialization
  63.   t = t0;
  64.   // the following step parameters are used in ODE45
  65.   // hmax = (tfinal - t)/5;
  66.   // hmin = (tfinal - t)/20000;
  67.   // h = (tfinal - t)/100;
  68.   // The following parameters were taken because the integrator has
  69.   // higher order than ODE45. This choice is somewhat subjective.
  70.   hmax = (tfinal - t)/2.5;
  71.   hmin = (tfinal - t)/10000;
  72.   h = (tfinal - t)/50;
  73.   y = y0[:];
  74.   f = y*zeros (1,13);
  75.   tout = t;
  76.   yout = y.';
  77.  
  78.   if (!exist (Uoutput)) 
  79.   {
  80.     output = 1;
  81.     lout = << [tout, yout] >>;
  82.   else
  83.     output = 0;
  84.     lout = << Uoutput (t, y).' >>;
  85.   }
  86.   I = 2;
  87.  
  88.   tau = tol * max( [norm(y, "I"), 1]);
  89.  
  90.   if (trace) { [t, y[:].'] ? }
  91.  
  92.   // The main loop
  93.  
  94.   while ((t < tfinal) && (h >= hmin))
  95.   {
  96.     if (t + h > tfinal) { h = tfinal - t; }
  97.  
  98.     // Compute the slopes
  99.     f[;1] = Ufun (t, y)[:];
  100.     for (j in 1:12)
  101.     {
  102.       f[;j+1] = Ufun (t+alpha[j]*h, y+h*f*beta[;j])[:];
  103.     }
  104.  
  105.     // Truncation error term
  106.     gamma1 = h*41/840*f*psi;
  107.  
  108.     // Estimate the error and the acceptable error
  109.     delta = norm (gamma1,"I");
  110.     tau = tol*max ([norm (y,"I"),1.0]);
  111.  
  112.     // Update the solution only if the error is acceptable
  113.     if (delta <= tau)
  114.     {
  115.       t = t + h;
  116.       y = y + h*f*chi;
  117.       if (output) 
  118.       {
  119.         lout.[I] = [t, y.'];
  120.       else
  121.         lout.[I] = Uoutput (t, y).';
  122.       }
  123.       I++;
  124.     }
  125.     if (trace) { [t, y[:].']; }
  126.  
  127.     // Update the step size
  128.     if (delta != 0.0)
  129.     {
  130.       h = min ([hmax, 0.8*h*(tau/delta).^pow]);
  131.     }
  132.   }
  133.  
  134.   if (t < tfinal)
  135.   {
  136.     fprintf ("stderr", "SINGULARITY LIKELY, at t = %d\n", t);
  137.   }
  138.  
  139.   return collapse (lout);
  140. };
  141.  
  142. //
  143. // Collapse a LIST of matrices into a single matrix.
  144. //
  145.  
  146. collapse = function ( LIST )
  147. {
  148.   local (i, m, row, col);
  149.  
  150.   row = size (LIST.[members (LIST)[1]])[1];
  151.   col = size (LIST.[members (LIST)[1]])[2];
  152.   m = zeros (length (LIST)*row, col);
  153.  
  154.   for (i in 1:length (LIST))
  155.   {
  156.     m[i;] = LIST.[i];
  157.   }
  158.   return m;
  159. };
  160.